Statistical Methods in Medical Research

Selection bias is a common concern in causal inference. Two mechanisms can give rise to selection bias: collider bias or stratifying on an effect modifier of the exposure-outcome causal effect. Research into the impact of selection bias on instrumental variable (IV) analyses has focused primarily on collider bias. Here, we show that IV estimates can be biased if selection is an effect modifier (or caused by an effect modifier) of the exposure-outcome effect. In a simulation study, we investigate the size and direction of the bias when estimating effects on the mean difference or odds ratio scale. Bias can also occur due to heterogeneity in the instrument-exposure relationship, although its magnitude is typically smaller. Finally, we discuss methods to deal with this bias in practice, including inverse probability weighting (IPW).

Purpose: Videourodynamics (VUDS) are used to assess bladder dysfunction and risk of upper urinary tract deterioration in patients with spina bifida (SB), but interpretation is subject to interrater variability. Machine learning models have previously classified bladder dysfunction severity and predicted incident hydronephrosis using VUDS data. This study evaluated the performance of adapted versions of these models in an independent external cohort. Materials and Methods: VUDS data were collected from SB patients who underwent VUDS at a single institution between 2016 and 2025. Previously developed machine learning frameworks were adapted for implementation in the independent dataset. Three models were utilized: a deep learning convolutional neural network model using pressure-volume data, a deep learning imaging model using fluoroscopic imaging data, and an ensemble model that averaged the risk data from the pressure-volume and fluoroscopic data. The models were evaluated for prediction of incident hydronephrosis and classification of bladder dysfunction severity compared to expert pediatric urologist reviewers. Results: 70 patients were included in the hydronephrosis cohort, of whom 13 (18%) developed incident hydronephrosis. The ensemble model outperformed individual modalities, with a concordance index of 0.77 (pressure-volume only: 0.70, imaging only: 0.72) and an overall AUROC of 0.75 (pressure-volume: 0.67, imaging: 0.74). For bladder dysfunction classification in 95 VUDS studies, the ensemble model achieved 71% accuracy compared to 60% for pressure-volume data alone and 66% for imaging alone, with no substantial disagreements with expert reviewers. Conclusions: The machine learning models that were developed to predict incident hydronephrosis and classify the severity of bladder dysfunction using VUDS data were able to be adapted to an external, independent dataset and perform with similar discrimination and accuracy as previously reported.

Estimating the correlation between subject-specific rates of change is a common objective in longitudinal studies. A natural approach is to fit ordinary least squares regressions per subject and correlate the resulting slope estimates directly. Under complete data, this slopes-first estimator is algebraically equivalent to a contrast-based plug-in estimator derived from the joint variance-covariance matrix of the two signals across time. Under independent dropout, however, the slopes-first approach implicitly applies heterogeneous weights across missing-data patterns, introducing bias whenever the cross-covariance between the two signals is not constant over time. We formalise this result and propose an unbiased alternative that estimates the joint covariance matrix before applying the contrast. We further extend our results to include covariates as well. Simulations confirm that the slopes-first estimator is biased under non-constant cross-covariance while the proposed alternative remains unbiased. Further, we show in a real-world example that the two estimators can lead to different scientific conclusions and recommend that researchers working with longitudinal data inspect the missing data pattern to make an informed decision about the estimator of choice.

Predicting implantation outcomes in vitro fertilization (IVF) is challenging because clinically relevant information becomes available at different points during treatment. In this study, we developed a stage-aware, leakage-conscious, and explainable machine-learning framework to predict implantation using a public IVF dataset. {beta}-HCG was excluded to avoid potential downstream target leakage. The data were divided into training, validation, and independent test sets. Several machine-learning models, including Logistic Regression, Random Forest (RF), Extra Trees, XGBoost, and Hist Gradient Boosting, were optimized using stratified cross-validation. A fixed RF model was then evaluated at three stages of the IVF process, using baseline couple characteristics, gamete and fertilization variables, and embryological variables, respectively. Model performance improved as additional treatment-related information became available, with the largest improvement observed after the inclusion of embryological variables. The final RF model showed good predictive performance on the independent test set. The number of embryos was consistently identified as the most important predictor across RF, permutation, and SHAP analyses. Overall, the findings suggest that IVF implantation prediction becomes more accurate as later-stage treatment information is incorporated, particularly information related to embryological characteristics. The proposed approach also provides interpretable predictions while reducing the risk of information leakage. Further evaluation using independent clinical datasets is needed before the framework can be considered for clinical use.

In systematic review, the measures reported by the primary studies may not always be the measures sought by the systematic review investigators. Motivated by an systematic review of interventions to reduce meat consumption, we discuss here methods to infer relative risks and confidence intervals from randomized controlled trials that reported differences, possibly in the presence of repeated measures or multiple interventions.

Background: Stepped-wedge cluster randomized trials (SW-CRTs) are increasingly used when interventions need to be implemented sequentially across clusters. Despite extensive developments in the design and analysis of superiority SW-CRTs, methods for planning noninferiority SW-CRTs remain limited. Methods: We propose swNIsim, a simulation-based framework for power calculation in noninferiority SW-CRTs using generalized linear mixed-effects models (GLMMs). The proposed approach accommodates both continuous and binary outcomes and supports practical design features, including unequal numbers of clusters across sequences, varying cluster-period sizes, delayed intervention effects, and analyses with or without adjustment for secular time effects. Noninferiority margins can be specified using several commonly applied strategies. Results: Through extensive simulation studies, we examined the impact of key design parameters on statistical power, including the number of periods, the number of clusters per sequence, the number of individuals per cluster-period, and the degree of between-cluster variability. For a fixed total sample size, designs with more time periods consistently achieved higher power. Increasing the number of clusters per sequence produced modest gains in power, whereas greater between-cluster variability reduced power. Type I error rates were well controlled across the scenarios examined. Balanced allocation of clusters across sequences was more efficient than unbalanced allocation schemes. Two real-world SW-CRT examples were used to illustrate application of the method. Conclusions: The proposed swNIsim provides a flexible and practical approach for power calculation in noninferiority SW-CRTs. By accommodating a wide range of design configurations and outcome types, and through implementation as a freely available R package, swNIsim may facilitate the planning and design of noninferiority stepped-wedge trials.

Existing prediction models for obstetric anal sphincter injuries (OASIS) often rely on traditional regression, lack external validation, and, when fitted on multicenter data, fail to account for hospital-level clustering. We aimed to develop and temporally validate mixed-effects machine learning models for OASIS and compare them with traditional regression models. A secondary aim was to evaluate whether macrosomia ([&ge;]4000 g) could substitute for continuous birth weight. We conducted a retrospective nationwide cohort study of singleton, term ([&ge;]37 weeks), cephalic vaginal deliveries in the Dutch Perinatal Registry. Data from 2016 to 2019 were used for model development with three-fold cross-validation; data from 2020 were reserved for temporal validation. Candidate predictors included 18 demographic, obstetric, and neonatal factors. We developed six models: mixed-effects gradient boosting (MEGB), mixed-effects LASSO (ME-LASSO), and backward stepwise logistic regression, each with either birth weight or macrosomia. Discrimination was assessed by the area under the receiver operating characteristic curve (AUC), calibration by calibration plots, and prediction accuracy by the Brier score. Among 614,255 deliveries, OASIS occurred in 2.6%. In temporal validation including macrosomia, MEGB achieved the highest discrimination (AUC 0.727; 95% CI 0.719-0.735), significantly outperforming ME-LASSO (AUC 0.719; 95% CI 0.710-0.727; P<.0001) and stepwise logistic regression (AUC 0.714; 95% CI 0.706-0.722; P<.0001). Including birth weight marginally improved performance, with macrosomia serving as an effective proxy. Consistently important predictors were nulliparity, prolonged second stage, mediolateral episiotomy, advanced maternal age, operative vaginal delivery, macrosomia, and previous cesarean delivery. Uncalibrated MEGB showed poorer calibration, but isotonic recalibration substantially improved it. After calibration, all models showed excellent agreement between predicted and observed risks. MEGB provided improved discrimination and comparable calibration relative to ME-LASSO and stepwise logistic regression. Replacing continuous birth weight with a binary macrosomia indicator resulted in minimal performance loss, suggesting it may suffice for clinical risk prediction across all three modeling approaches.

Multiview data analysis has become increasingly important for integrating heterogeneous information from multiple data sources. Cooperative learning provides a flexible framework for multiview integration by bridging early and late fusion strategies. However, its use of squared loss makes it sensitive to outliers. In this study, we propose robust cooperative learning based on the soft-IPOD framework. The proposed method introduces an observation-specific outlier absorption vector into the cooperative learning objective, thereby reducing the influence of outlying observations while preserving information sharing across views. The resulting optimization problem can be formulated as a weighted Lasso-type problem using an augmented design matrix and can be naturally extended to more than two data views. Numerical simulations and a real-data application to labor onset prediction demonstrate the predictive performance of the proposed method in the presence of outliers.

Host immunity acquired through prior infections is a key factor shaping the spread of emerging variants. However, accounting for complete infection histories in mathematical models is challenging because the number of possible histories grows exponentially with the number of variants. Here, we present a flexible framework for modeling the co-circulation of multiple variants emerging at arbitrary times while incorporating cross-immunity and variant-specific characteristics. We introduce three tractable approximations: truncation, which retains only infection histories up to a given number of strains; one variant at a time, which is exact when variants circulation does not overlap; and two variants at a time, which is exact when at most two strains co-circulate. We systematically evaluate these approximations across a range of cross-immunity levels and inter-variant emergence delays. Truncation approximations perform best under short emergence delays and strong cross-immunity, whereas the other approaches outperform them at longer delays. When applied to SARS-CoV-2 variant dynamics, the second-order truncation and the two variants at a time approximation accurately capture incidence patterns, whereas the other approaches perform well only for specific subsets of variants. These results identify the parameter regimes in which simplifying assumptions remain valid, provide guidance for model selection, and highlight the key role of infection history in shaping epidemic trajectories.

Simulation models are widely used to support counterfactual inference in settings where direct experimental validation is difficult or infeasible. However, for phenomenological models calibrated to observational time series, goodness-of-fit does not necessarily establish that the structural assumptions underlying such inferences are uniquely supported by the data. Here, we examine this issue using a published COVID-19 transmission model for Japan that incorporates empirical inputs including human mobility, variant prevalence, and vaccination and has been used to estimate the effects of alternative vaccination scenarios. Using the publicly archived code together with a reconstructed substitute for a restricted-access input file required by the released implementation, we first reproduced the published model behaviour. We then performed two complementary computational stress tests. First, the vaccination-related input was suppressed and the remaining parameters were re-estimated against the original epidemic trajectory. Second, as a complementary test of trajectory-level flexibility, the model was refitted to synthetic target series with qualitatively different temporal structures while retaining the empirical mobility and variant inputs. The model retained substantial fitting capacity in both tests: suppression of the vaccination-related input caused only a limited deterioration in fit to the original trajectory, while qualitatively different synthetic trajectories could also be reproduced through parameter adjustment. Variation in fitted parameters across these conditions was also consistent with compensatory parameter behaviour. These results do not imply that vaccination had no effect or that the original model is invalid. Rather, they provide empirical evidence that, in this case, goodness-of-fit to the observed epidemic trajectory alone does not uniquely support the structural interpretation required for vaccination-related counterfactual inference. The study illustrates the value of reproducible computational stress testing as a complement to, rather than a substitute for, formal identifiability or model-comparison analyses.

Mechanistic models of infectious disease transmission provide a framework for understanding and predicting disease dynamics, but statistical inference is often computationally prohibitive. This challenge is especially pronounced for models with high-dimensional structures and many unknown parameters, where trajectory matching may require thousands of evaluations of an expensive likelihood function. We developed a surrogate-assisted optimization framework that uses Bayesian optimization to reduce the number of evaluations of the underlying mechanistic model required for parameter inference. We evaluated the approach using simulated fitting problems based on age-structured susceptible-exposed-infected-recovered models with 4, 8, and 16 age classes and a 210-compartment age-structured rotavirus transmission model. Applying surrogate-assisted optimization substantially reduced the number of likelihood evaluations required to identify good parameter estimates, with larger improvements as model dimensionality increased. In the rotavirus model, the surrogate Bayesian optimization approach reached the benchmark after approximately 0.7 hours, compared with approximately 9.4 hours using traditional Nelder-Mead optimization. These results demonstrate that surrogate-assisted optimization can substantially reduce the computational burden of trajectory matching for complex mechanistic transmission models, particularly for higher-dimensional inference problems.

Non-exponential residence times give rise to history-dependent or memory effects in epidemic dynamics, yet their mathematical formulations remain fragmented. Here we develop a unified framework for a history-dependent SEIR model with general latent- and infectious-period distributions. We first establish the equivalence between the stage-age and renewal formulations and show that both admit a unified transition-flow representation. We then introduce the concept of a transition-flow multiplier, which provides a direct link between residence-time distributions and transition flows. When the multiplier is a weighted sum of shifted power functions, the model admits an exact representation in terms of Riemann-Liouville or tempered Riemann-Liouville fractional operators, providing a first-principles foundation for fractional epidemic models. Finally, using Euler-Lotka equations and stochastic-order arguments, we characterize how the latent- and infectious-period distributions shape epidemic growth and decay. A longer latent period slows both early growth and late decay, while, at fixed means, greater variability in the latent period accelerates early growth but slows late decay, whereas greater variability in the infectious period slows both.

Background and Objective: Adding belzutifan to adjuvant pembrolizumab improved disease-free survival (DFS) in LITESPARK-022, but its durability beyond the 1-year treatment is unknown. Methods: Following a registered plan, we reconstructed individual patient data from published Kaplan-Meier curves of ten phase 3 adjuvant renal cell carcinoma trials. Within each trial, we estimated restricted mean survival time differences ({Delta}RMST), hazard ratios (HRs) before and after planned treatment ended, and DFS differences. Key Findings and Limitations: In LITESPARK-022 (central review), the 24-month {Delta}RMST was 1.34 months (95% confidence interval 0.73-1.95), and the belzutifan arm remained ahead after month 12 ({Delta}RMST increment 0.74 months/yr, 0.32-1.19). The HR was 0.55 (0.43-0.71) before month 12 and 1.16 (0.87-1.54) afterwards (post/pre ratio 2.10); investigator assessment showed similar attenuation. Similar ratios occurred in KEYNOTE-564 at its first interim analysis (2.02), whose benefit persisted, and in SORCE and PROTECT (2.09, 2.59), whose benefit faded. From months 12 to 30 (post hoc), the DFS difference changed by -3.7 percentage points (-8.1 to 0.8) with central review and 0.6 (-4.0 to 5.0) with investigator assessment. Data were reconstructed, and follow-up was short. Conclusions and Clinical Implications: Relative hazard reductions with belzutifan were greater before month 12; durability of the absolute DFS advantage is uncertain and assessment-dependent, so single-HR projections should be avoided.

Accurate prediction of individual medical outcomes is essential for optimizing treatment allocation amid rising costs, coverage denials, and limited clinical resources. Traditional predictive models, including linear regression and neural networks, rely on average effects and cannot tailor predictions to the specific circumstances of individual cases. We present relevance-based prediction (RBP), a model-free and task-specific method that predicts outcomes as weighted averages of observed cases, with weights determined by a rigorously defined and theoretically justified measure of relevance. Unlike model-based methods that rely on fixed calibrated parameters, relevance-based prediction revisits the original data for each prediction and customizes both the cases and variables used. We illustrate RBP by applying it to predict opioid treatment outcomes and demonstrate that it provides case-specific insights unavailable from conventional models, including how each prior case informs a prediction, how each variable affects each predictions reliability and value, and how reliable each prediction is before it is made. These individualized insights may prevent misleading average-based decisions and reduce harmful or suboptimal treatment.

Doubly robust survival estimators require nuisance models that the available information can support. We examine design and estimator dependence in TCGA breast cancer data. Treatment indicators disagree for 337 of 1076 patients; 67.4% have unobserved five-year binary outcomes; and receptor imputation expands the triple-negative subgroup from 156 directly measured patients to 226. We compare chemotherapy initiation by day 90 with non-initiation, using adjusted restricted mean survival time (RMST) over the next 1095 days. Weighting-only estimates are 55 days (95% bootstrap interval: -24, 133) in triple- negative disease (16 deaths) and 38 days (12, 72) in the whole cohort (72 deaths). A Cox outcome-model benchmark yields 108 and 35 days, respectively, demonstrating material model dependence in the smaller population. A parameter-to-information ratio summarizes the burden on complete-outcome regressions within cross-fitting. Its simulation-calibrated warning level is 0.31 (0.28 - 0.35) for the original regime; additional simulations examine covariate-dependent censoring, delayed benefit and nonlinear misspecification. The ratio does not certify estimator validity or replace weight diagnostics. The contribution is a reproducible workflow that separates design, nuisance-model support and estimator stability before interpreting observational survival contrasts.

Background. Classical macroparasite models commonly represent worm burden heterogeneity through prescribed Poisson or negative-binomial (NB) distributions, while transmission dynamics are formulated in terms of aggregate quantities such as mean worm burden and prevalence. We revisit these assumptions using a stratified worm burden (SWB) framework that explicitly tracks the distribution of parasite burdens across heterogeneous host classes and allows demographic turnover, transmission, and interventions to be analyzed within a unified model. Methods. We derive equilibrium distributions of the SWB system under homogeneous and heterogeneous infection processes and compare them with the classical Poisson and NB distributions used in macroparasite epidemiology and control policymaking. The framework is then extended to include host demography, coupled human-snail transmission dynamics, and periodic mass drug administration (MDA). Corresponding SWB and NB-type MacDonald models are calibrated to comparable endemic prevalence and mean worm burden, allowing direct comparison of potential transmission breakpoints and intervention outcomes. Results. In the absence of demographic turnover, the SWB framework naturally generates Poisson and NB equilibrium distributions as limiting cases of homogeneous and heterogeneous infection processes. Incorporating demographic recruitment and mortality alters these benchmark distributions, producing endemic burden structures that differ systematically from their classical counterparts. Analysis of coupled human-snail systems reveals that transmission breakpoints and bistable equilibria occur only within a relatively restricted region of SWB parameter space. In contrast, comparable NB-MacDonald formulations exhibit breakpoint behavior over substantially broader parameter ranges with the net effect of overestimating the likelihood of interruption of transmission. Although calibrated models can reproduce near identical endemic prevalence and mean worm burden, they often generate markedly different responses to MDA, including differences in rebound dynamics, persistence, and elimination thresholds. Conclusions. The SWB framework links classical statistical descriptions of worm aggregation to mechanistic burden-resolved population dynamics, clarifying how endemic distributions emerge and how they are modified by demography and control interventions. Our results suggest that transmission breakpoints are not a universal feature of helminth systems but depend strongly on model structure, and that predictions of elimination and post-treatment recovery may be considerably more model-dependent than implied by prevalence-based calibration alone. The SWB approach therefore provides a flexible alternative framework for evaluating parasite transmission dynamics and predicting the outcomes of control strategies.

Importance. Observational studies dominate evidence in cardiovascular outcomes research, yet no widely adopted reporting standard requires authors to disclose how much of an observed association is attributable to systematic error. The negative-control calibration framework of the Observational Health Data Sciences and Informatics consortium returns a binary answer: a calibrated P-value either crosses 0.05 or it does not. That answer conflates "is there a signal" with "is the signal credible." Objective. To define the Bias Attribution Fraction (BAF), a continuous estimator in [0, 1] that quantifies the fraction of an observed log relative risk that is attributable to systematic bias, together with its 95% credible interval, and to evaluate whether reporting calibrated P-value plus BAF plus interval width produces more defensible decisions than reporting the calibrated P-value alone. Design. Three-part study: (i) metric definition and interval construction; (ii) Monte Carlo simulation across 1920 conditions by 1000 repetitions = 1,920,000 simulated target estimates (a primary grid plus a sign-flipped companion grid of identical size, pooled for all analyses), with truth known by construction; (iii) two-arm beta-blocker target trial emulation in MIMIC IV, adhering to the TARGET guideline, with targeted maximum likelihood estimation (TMLE) and 12 negative-control outcomes. Setting. MIMIC IV (critical care) and MIMIC IV ECG, 2008 to 2019. Participants. Adults with a first intensive care unit stay and a recorded 12-lead ECG. Exposures. Two clinical questions are evaluated: (1) reaching at least 50% of guideline-recommended target dose of an index beta-blocker within 7 days of admission, sustained to 12 months, versus lower or no exposure; (2) receiving guideline-directed medical therapy (GDMT) with at least two of three foundational drug classes at adequate dose for 12 months versus fewer. Main Outcomes and Measures. One-year all-cause mortality. Primary analytic outputs: (i) the calibrated P-value from a negative-control empirical null distribution with 12 outcomes; (ii) BAF-hat with a Markov chain Monte Carlo 95% credible interval; (iii) a verdict drawn from a truth-calibrated lookup table mapping the BAF-hat by interval-width cell to a probability of being bias-dominated. Results. Across 1,920,000 simulated repetitions, the calibrated P-value alone yielded a misuse rate (declared usable while truth was bias-dominated) of 36.02%. Adding the BAF point estimate reduced misuse to 16.48%, and requiring the whole credible interval of BAF to lie below 0.5, clear of the bias-dominated zone, reduced misuse to 10.16%. Requiring in addition a calibrated probability of bias dominance below 15% reduced misuse to 3.02% at the cost of declaring 22.1% of estimates usable (versus 61.8% for the calibrated P-value alone). Read from the bias-dominated view, meaning that withholding a repetition is a positive call and being truly bias-dominated is a positive truth, the BAF estimator achieved an area under the receiver operating characteristic curve of 0.922, compared with 0.787 for the calibrated P-value. In MIMIC IV, question 1 produced a calibrated P-value of 0.894 and BAF-hat = 0.875 (95% credible interval 0.385 to 0.992); the apparent one-year mortality reduction of 16.9% was entirely attributable to systematic bias and the association was not usable as effect evidence. Question 2 also failed layer 1, with a calibrated P-value of 0.1026, and BAF-hat = 0.460 placed it in the "competitive, no verdict" zone. Conclusions and Relevance. Pairing the calibrated P-value with BAF and its 95% credible interval turns a binary credibility judgement into a three-channel report (signal existence, effect size, and bias fraction). Despite TARGET-adherent design, doubly robust TMLE, and 12 negative controls, both clinical questions in this study still failed to qualify as effect evidence, with question 1 demonstrating complete bias domination and BAF quantifying the systematic error that remained. The reporting standard proposed here can be applied to any observational study that includes negative controls, without requiring additional data.

Background: Excess mortality estimates depend on an expected baseline, an unobservable quantity typically derived from models with prescribed functional forms and short reference periods that exclude exceptional events. Singular Spectrum Analysis (SSA) offers a nonparametric alternative but requires subjective component selection and is sensitive to outliers. We combined SSA with automated component selection and outlier interpolation to estimate long-term baselines and compile a catalog of excess mortality events (EMEs). Methods: We analyzed weekly all-cause mortality (2015-2024) for six metropolitan areas with contrasting climatic, demographic, and epidemiological contexts: Brussels, Antwerp, and Liege (Belgium), and Athens, Thessaloniki, and Larisa (Greece). SSA was applied to standardized log-mortality rates; components were selected using Monte Carlo SSA, and COVID-19 periods were interpolated iteratively. EMEs were defined as weeks exceeding a quasi-Poisson two-z-score threshold, then screened for coincidence with heat waves, cold waves, COVID-19 periods, and influenza epidemics. Results: Baselines were consistent within countries: annual cycles in Belgium, with an additional significant semiannual cycle in Greece reflecting summer mortality. We cataloged 170 EMEs spanning 456 weeks and 14,210 excess deaths. Pandemic weeks accounted for 77% of excess mortality, and 96 weeks coincided with temperature extremes; 56 EMEs, representing 3.7% of excess mortality, remained unmatched. Conclusions: SSA yields reproducible, low-assumption baselines that capture country-specific seasonality. Because flexible baselines can absorb recurrent seasonal mortality, component selection is a substantive choice that should be reported explicitly, together with its sensitivity.

Understanding the relationship between weather conditions and infectious disease transmission has important public health implications; nevertheless, the functional form of this relationship often remains poorly understood. Humidity has previously been suggested to have a U-shaped effect on the transmission of respiratory pathogens, with optimal transmission occurring at both low and high levels. Here, to assess the impact of this functional form and of its potential misspecification, we first simulate transmission rates based on the U-shaped relationship and then test whether an exponential function accurately estimates these rates across a range of climates. We further use a transmission model to simulate case count data for a range of hypothetical pathogens, comparing outbreak dynamics under the U-shaped forcing function to those under the best-fit exponential approximation. In our simulations, the model fit for the exponential function is climate-dependent, performing best in very dry and very humid environments. In locations where specific humidity values centre around the vertex of the U-shaped function, the exponential model underestimates the amplitude and frequency of incidence peaks. These findings demonstrate the potential consequences of assuming an incorrect functional form for weather-infectious disease relationships and highlight the need for more careful consideration of this functional form in modelling practice.

Purpose: Digital-twin frameworks can potentially offer a pathway toward individualized treatment but remain largely untested in radiation oncology. We present a causal digital-twin model for head-and-neck (H&N) cancer that estimates patient-specific survival benefit from chemoradiotherapy (CRT) versus radiotherapy (RT) alone, stratified by HPV status. Methods: Using the public RADCURE cohort (n = 3,346), we developed a causal digital-twin framework in which patients were divided into HPV-positive and HPV-negative groups and assigned to training (RT start year [&le;] 2007) and test sets (RT start year > 2007). Bayesian-network structure learning identified each cohort's Markov Blanket around the 2-year survival endpoint (SVy2), yielding the minimal causally relevant feature set. A causal-forest model (double machine learning, DML framework) estimated the individualized treatment effect {tau}=P(SVy_2|CRT,X)-P(SVy_2|RT,X), with {tau} denoting the model-estimated effect. Overlap restrictions were applied to mitigate extrapolation. Predictive discrimination (AUC, C-index), average treatment effects (ATE), individualized effects (ITE), and decision-curve analyses were computed. Results: For HPV-negative disease, internal evaluation demonstrated strong discriminative performance (test AUC = 0.72; C-index = 0.81, p < 0.01). The estimated test-set ATE favored CRT with a modest but positive improvement in SVy2 (mean {tau} = 0.13). Patterns in {tau} were clinically coherent, including increased benefit in node-positive and good-performance patients and reduced benefit in high-burden or poor-performance presentations. Decision-curve analysis showed that a {tau} -guided selective-CRT strategy provided higher net benefit than uniform treatment policies. In HPV-positive disease, treatment-effect heterogeneity was attenuated, consistent with its more uniformly favorable prognosis. Conclusion: This study demonstrates a proof-of-concept causal digital-twin framework for interpretable, patient-specific survival simulation in H&N cancer. The model captures clinically plausible treatment-effect heterogeneity and establishes a foundation for multi-cohort validation and future integration into individualized treatment decision support.